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Method  of  Lines  Transpose  an  Implicit  Vlasov  Maxwell  Solver  for  Plasmas 


ANDREW  CHRISTLIEB 

Abstract.  We  propose  a  new  particle-in-cell  (PIC)  method  for  the  simulation  of  plasmas 
based  on  a  recently  developed,  unconditionally  stable  solver  for  the  wave  equation.  This 
method  is  not  subject  to  a  CEL  restriction,  limiting  the  ratio  of  the  time  step  size  to  the 
spatial  step  size,  typical  of  explicit  methods,  while  maintaining  computational  cost  and  code 
complexity  comparable  to  such  explicit  schemes.  We  describe  the  implementation  in  one 
and  two  dimensions  for  both  electrostatic  and  electromagnetic  cases,  and  present  the  results 
of  several  standard  test  problems,  showing  good  agreement  with  theory  with  time  step  sizes 
much  larger  than  allowed  by  typical  CEL  restrictions. 


1.  Introduction 

This  is  a  final  report  for  the  AFOSR  grant  “Method  of  Lines  Transpose  an  Implicit  Vlasov 
Maxwell  Solver  for  Plasmas”.  The  contact  number  was  FA9550-11-1-0281.  The  grant  was 
concerned  with  the  development  of  a  new  approach  for  developing  A-sable  high  order  waves 
solvers  for  combination  with  particle  methods.  The  work  on  the  wave  solver  has  been  pub¬ 
lished  and  can  be  found  in  [8,  7,  9].  In  this  report,  we  detail  our  first  attempt  at  turning  this 
wave  solver  into  an  implicit  method  for  particle-in-cell,  as  well  as  go  over  a  host  of  example 
problems  in  ID  and  2D.  This  work  was  in  conjunction  with  two  former  post  docs,  funded 
by  this  grant.  Dr.  Yaman  Giiglu  and  Dr.  Matt  Caulesy,  as  well  as  my  student,  Mr.  Eric 
Wolf,  who’s  thesis  is  the  subject  of  this  work.  In  addition,  we  have  ben  collaborating  with 
Dr.  Matthew  Bettencourt  from  Sandia  National  Lab,  who  has  ben  part  of  the  large  scale 
effort  at  Sandia  to  develop  an  implicit  particle  method.  These  ideas  can  be  extended  to 
methods  for  Pseudo  Differential  Operators  (PDO),  and  my  current  student,  Ms.  Hana  Cho, 
who  is  starting  to  work  in  this  area. 

Collisionless  plasmas  -  systems  of  charged  particles  interacting  through  electromagnetic 
fields  -  are  modelled  by  the  Vlasov-Maxwell  system  of  partial  differential  equations  (PDEs), 
which  couple  Maxwell’s  equations,  describing  the  evolution  of  the  electric  and  magnetic 
fields  E  and  B,  to  Vlasov  equations,  a  type  of  hyperbolic  PDE  describing  the  evolution  of 
the  phase-space  distribution  functions  (DFs)  fs  of  the  various  species  s  of  charged  particles. 
Particle-in-cell  (PIC)  methods  [4,  18,  38],  in  development  and  use  since  the  1960s  and  a 
primary  tool  in  the  computer  simulation  of  plasmas,  combine  an  Eulerian  description  of  the 
fields  with  a  Lagrangian  description  of  the  DFs;  that  is,  fields  are  evolved  on  a  fixed  mesh, 
while  DFs  are  represented  by  moving  particles  whose  trajectories  are  characteristics  of  the 
corresponding  Vlasov  equation.  Thus,  PIC  methods  require  a  method  to  compute  the  fields 
on  a  mesh  and  a  method  to  compute  particle  trajectories,  as  well  as  interpolation  tools  to 
provide  for  their  coupling.  This  work  focuses  on  a  new  method  for  the  computation  of  the 
fields,  along  with  associated  interpolation  techniques. 
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Under  the  Lorenz  gauge  condition,  Maxwell’s  equations  reduce  to  uncoupled  wave  equa¬ 
tions  for  the  scalar  and  vector  potentials,  $  and  A.  Recently,  a  novel  method  for  the  solution 
of  the  wave  equation  has  been  developed  [8,  7,  9],  based  on  the  Method  of  Lines  Transpose 
(MOLT),  dimensional  splitting  and  an  efficient  ID  integral  solution  method,  which  is  un¬ 
conditionally  stable  (or  A-stable)  -  that  is,  it  is  not  subject  to  the  Courant-Friedrichs-Lewy 
(CFL)  restriction  limiting  the  ratio  of  the  temporal  step  size  to  the  spatial  step  size,  typical 
of  widely  used  explicit  methods.  In  this  work,  we  apply  this  method  to  the  uncoupled  wave 
equations  for  $  and  A  to  solve  Maxwell’s  equations  with  a  method  comparable  in  computa¬ 
tional  cost  and  complexity  of  code  to  explicit  methods  such  as  the  well-known  Yee  scheme 
[39,  36],  but  without  introducing  a  CFL  restriction  based  on  the  speed  of  light  as  in  such 
explicit  methods.  In  conjunction  with  an  appropriate  description  of  particles,  we  seek  to 
develop  a  PIC  method  that  retains  the  simplicity  of  explicit  finite-difference-based  methods 
while  eliminating  this  CFL  restriction. 

We  demonstrate  the  application  of  our  method  to  both  electrostatic  and  electromagnetic 
problems.  In  the  non-relativistic,  zero-magnetic  field  limit,  it  is  typical  to  make  the  elec¬ 
trostatic  approximation,  E  =  —  V<F,  —  =  p/eo,  simplifying  the  Vlasov-Maxwell  system 

to  the  Vlasov-Poisson  system.  Correspondingly,  in  this  work  we  consider  the  same  non- 
relativistic,  zero-magnetic  field  limit,  and  drop  A  and  the  corresponding  wave  equations 
from  our  model.  When  particle  velocities  are  small  compared  to  the  speed  of  light,  we  ar¬ 
gue  that  this  model,  which  we  term  the  Vlasov- Wave  model  due  to  the  replacement  of  the 
Poisson  equation  with  a  wave  equation,  will  agree  approximately  with  the  usual  electrostatic 
Vlasov-Poisson  model.  We  present  numerical  results  applying  our  method  to  several  stan¬ 
dard  electrostatic  test  problems  in  one  and  two  dimensions,  showing  agreement  with  the 
predictions  of  linear  theory  for  the  electrostatic  model.  We  apply  our  method  to  a  fully 
electromagnetic  beam  pinch  problem  in  two  dimensions.  It  is  well  known  that  an  electro¬ 
magnetic  PIC  method  must  satisfy  a  discrete  form  of  Gauss’  law  (V  ■  E  =  p/eo)  to  prevent 
serious  numerical  errors  related  to  the  violation  of  charge  conservation  [26] .  In  this  work,  we 
obtain  solutions  that  satisfy  exactly  discrete  forms  of  Gauss’  law  for  the  electric  field  and  the 
divergence-free  condition  for  the  magnetic  field  through  a  staggered  grid  approach,  adapted 
from  the  well-known  Yee  grid  [39],  with  a  Poisson  equation  formulation  for  the  scalar  po¬ 
tential.  In  addition  to  eliminating  the  GFL  restriction,  the  wave  solver  method  used  in  this 
work  offers  the  handling  of  complex  boundary  geometry  in  a  Gartesian  grid  without  using  a 
staircasing  approximation  [7,  37],  and  can  be  extended  to  higher-order  accuracy  [9],  features 
that  will  be  incorporated  into  our  PIG  method  in  future  work. 

Building  on  prior  work  on  unconditionally  stable  ADI-FDTD  schemes  for  Maxwell’s  equa¬ 
tions  [41,  40,  28,  23],  unconditionally  stable  ADI-FDTD  methods  for  Maxwell’s  equations 
that  conserve  the  divergence  of  the  discrete  electric  and  magnetic  fields  were  incorporated 
into  PIG  methods  in  [35].  Future  work  will  compare  our  method  to  this  approach  in  the 
context  of  PIG  methods. 

In  most  PIG  methods,  further  stability  restrictions  also  apply,  the  most  restrictive  often 
being  the  need  to  resolve  the  electron  plasma  period.  Our  method,  making  use  of  explicit 
algorithms  for  advancing  particles,  is  subject  to  such  restriction.  In  problems  where  time 
scales  much  longer  than  the  electron  plasma  period  are  of  primary  interest  and  dynamics  on 
the  scale  of  the  electron  plasma  period  can  be  safely  underresolved  -  for  instance,  in  certain 
problems  in  the  study  of  ion  dynamics  -  it  is  desirable  to  take  a  much  larger  time  step  than 
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prescribed  by  the  plasma  period  stability  restriction.  Since  years  ago  [14,  22],  and  with  a 
recent  resurgence  [11],  it  has  been  sought  to  develop  implicit  PIC  algorithms  that  are  not 
subject  to  the  stability  restriction  based  on  the  plasma  period.  An  ultimate  challenge  for  our 
method  would  be  synthesis  with  a  suitable  implicit  particle  integration  scheme  to  achieve 
a  practical  fully  implicit  method,  eliminating  the  stability  restriction  based  on  the  plasma 
period  as  well  as  the  field-based  CFL  restriction.  However,  that  is  beyond  the  scope  of  the 
present  work. 


2.  Physical  Models 

The  self-consistent  evolution  of  an  single  species  plasma  is  described  by  the  Vlasov-Maxwell 
system,  given  in  the  SI  system  of  units  by 


^  +  V  •  Vx/  +  -  (E  +  V  X  B)  ■  Vv/  =  0 
ot  m 


—  V  X  B  —  /ioJ 
V  ■  E  =  p/^o 


dB 

V  B 


-V  X  E 
0 


p(x,  t)=q  /(x,  V,  t)d\,  3  =  q  v/(x,  v,  t)  dv 


where  /(x,  v,  t)  is  the  phase-space  distribution  function,  q  is  the  charge  and  m  is  the  mass  of  a 
particle,  E(x,  t)  is  the  electric  field,  B(x,  t)  is  the  magnetic  field,  p(x,  t)  is  the  charge  density 
(and  pb  represents  a  static  uniform  background  charge  distribution),  J(x, f)  is  the  current 
density,  Eq  is  the  electric  permittivity  of  the  vacuum,  and  po  is  the  magnetic  permability 
of  the  vacuum.  (Boldface  variables  are  to  stand  for  vector  quantities,  while  nonboldface 
variables  are  to  stand  for  scalar  quantities.) 

In  the  limit,  in  an  appropriate  sense,  of  |v|/c  — ?■  0  (where  c  =  1/ ^EqPq  is  the  speed  of 
light  in  vacuum)  and  in  the  absence  of  magnetic  fields,  the  Vlasov-Maxwell  system  reduces 
to  the  Vlasov-Poisson  system: 

+  v-Vx/--V$-Vv/  =  0 
ot  m 


-A$  =  p/eq, 


p(x,f)  =  q  /(x,v,f)dv 


The  Vlasov-Maxwell  and  Vlasov-Poisson  models  have  been  widely  studied,  and  many 
numerical  methods  have  been  developed  to  solve  them,  such  as  electrostatic  and  electromag¬ 
netic  PIC  methods  [4,  18,  13,  27,  17],  as  well  as  Eulerian  methods  [12,  16,  31,  1].  Instead  of 
solving  either  of  these  systems  directly,  we  seek  to  develop  a  semi-implicit  approach  for  the 
fields  based  on  a  vector  potential  formulation.  Under  the  Lorenz  gauge  condition. 


1  ad) 

dt 


+  V- A 
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Maxwell’s  equations  reduce  to  uncoupled  wave  equations; 


1  ^2$ 

dt'^ 
1  d^A 


■  =  p/eo 

AA  =  /ioJ 


with  E  =  -V<h  -  ^  and  B  =  V  X  A  [21], 

In  the  electromagnetic  case,  we  will  solve  these  wave  equations  coupled  into  a  particle- 
in-cell  method  to  constitute  a  solution  of  the  usual  Vlasov-Maxwell  system.  In  the  non- 
relativistic,  zero-magnetic  field  limit,  we  drop  A  and  the  corresponding  wave  equations  and 
consider  a  quasi-electrostatic  Vlasov-Wave  model,  which  will  agree  with  the  electrostatic 
Vlasov-Poisson  model  when  |v|/c  <<  1.  This  condition  frequently  can  be  interpreted  as 
UpL  <<  c  for  a  relevant  physical  length  scale  L,  where  Up  is  the  (electron)  plasma  frequency. 
The  resulting  Vlasov- Wave  system  is  as  follows: 


^  +  v-Vx/--Vd>-Vv/  =  0 
ot  m 

/(x,v,f)civ 

The  purpose  of  this  work  is  the  numerical  solution  of  (nondimensionalized  forms  of)  these 
systems  with  PIC  methods  that  avoid  imposing  a  CFL  stablity  restriction  related  to  the 
speed  of  light  due  to  the  presence  of  the  wave  equations  for  $  and  A. 


1  ^2$ 


-A^  =  p/eo,  p(x,f)  =  g 


3.  Wave  Solver  Description 

3.1.  Origin  and  Mitigation  of  the  CFL  Restriction.  Consider  the  Cauchy  problem  for 
the  wave  equation: 


Utt  =  (?Au 

(1) 

M(a;,0)  =  g{x) 

(2) 

Ut{x,  0)  =  h{x). 

(3) 

To  see  the  origin  of  the  CFL  restriction  in  a  typical  explicit  finite  difference  method, 
consider  the  semi-discrete  scheme  obtained  by  substituting  the  centered,  second-order  finite 
difference  discretization  of  Uu  into  the  wave  equation.  We  obtain 


u{x,  t  +  At)  —  2u{x,  t)  +  u{x,  t  —  At) 
At2 


c^A‘^u{x,  t). 


Upon  Fourier  transforming  in  the  spatial  variable,  we  obtain 


t  +  At)  —  2u{^,  t)  -|-  t  —  At)  =  —c^At‘^\^\‘^u{^,  t)  (4) 

where  ^  is  the  spatial  frequency.  The  problem  is  that  a  high-frequency  perturbation  of  u 
-  such  as  that  introduced  by  truncation  error  -  will  be  amplified  by  the  symbol  — c^At^|^p, 
resulting  in  instability.  In  a  fully  discrete  method,  the  spatial  frequencies  are  bounded  by 
2/ Ax  (the  Nyquist  frequency),  meaning  that  cAt/Ax  must  be  sufficiently  small  to  prevent 
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amplification  of  high  frequencies.  In  terms  of  semi-discrete  Von  Neumann  analysis,  consider 
substituting  into  4  a  semi-discrete  solution  of  the  form  u{^,nAt)  =  A”  for  some  A  G  C.  We 
obtain  a  quadratic  equation  A^  —  (2  —  ^^)A  -|- 1  =  0  and  the  stability  condition  | A|  <  1  for  any 
root  A  of  this  equation,  where  z  =  cAt|^|.  It  is  easily  verihed  from  the  quadratic  formula 
that  the  roots  Aj  satisfy  |Aj|  <  1  for  i  =  1,2  if  and  only  if  \z\  <  2,  that  is  if  and  only  if 
cAt  <  ‘2/\^\.  For  the  ID  case,  inserting  the  Nyquist  frequency  2/ Ax  into  the  upper  bound 
gives  the  usual  CFL  stability  restriction  cAt  <  Ax. 

To  avoid  the  CFL  restriction  while  maintaining  an  explicit  update  equation,  we  modify 
the  symbol  of  the  Laplacian  to  bound  it  and  prevent  the  amplihcation  of  high  frequencies. 
Consider  the  one-dimensional  case.  We  look  for  bounded  rational  approximations  of  near 
.^  =  0.  One  option  is 


4"^  -I- 

We  can  write  this  modihed  semi-discrete  scheme  in  Fourier  space  as 

Q,2f2 

m(4,  t  +  At)-  2h(4,  t)  +  m(4,  t-  At)  =  -c^At^— — t) 

4  + 

Choosing  a  =  jS/cAt  for  some  properly  chosen  /?,  we  can  see  that  the  modihed  symbol 
is  uniformly  bounded  by  /3‘^  for  all  At,  so  that  choosing  /?  <  2  will  ensure  stability,  inde¬ 
pendently  of  At,  and  it  turns  out  that  /3  =  2  gives  optimal  accuracy  for  stable  schemes  of 
this  family  [7,  9].  Multiplication  in  Fourier  space  is  equivalent  to  convolution  in  physical 
space,  so  transforming  back  into  physical  space  gives  a  method  based  on  convolution  with 
the  inverse  Fourier  transform  of  the  modihed  symbol; 

u{x,  t  +  At)—2u{x,  t)  +  u{x,  t  —  At)  =  —l3‘^D[u]{x,  t) 

where  D[f]{x)  =  f  (d(x  —  x')  —  dx'  =  f{x)  —  f  ■ 

As  a  general  comment,  this  modihcation  of  the  symbol  of  the  Laplacian  introduces  error 
at  high  spatial  frequencies,  so  that  this  technique  may  not  be  suitable  for  problems  in  which 
extremely  accurate  propagation  of  waves  at  high  spatial  frequencies  (high  wavenumber)  is 
essential.  However,  in  many  problems  of  plasma  physics,  the  physics  is  dominated  by  ehects 
at  low  spatial  frequencies,  and  in  these  problems  this  technique  can  hnd  useful  application. 

A  detailed  discussion  of  numerical  methods  for  the  wave  equation  derived  from  the  per¬ 
spective  of  Fourier  multipliers,  including  higher-order  methods  and  the  extension  to  multiple 
dimensions,  can  be  found  in  [9].  In  the  following  section,  we  give  the  derivation  of  two  useful 
second-order  numerical  methods  through  alternative  means,  the  Method  of  Lines  Transpose, 
which  also  facilitates  the  extension  of  the  methods  to  bounded  domains  and  to  multiple 
dimensions. 

3.2.  Method  of  Lines  Transpose:  Derivation  of  Semi-Discrete  Schemes.  In  the 

Method  of  Lines  Transpose  for  the  solution  of  time- dependent  PDFs,  also  known  as  Rothe’s 
method  [33],  hnite  difference  discretizations  of  time  derivatives  are  substituted  into  the  PDE, 
resulting  in  a  boundary  value  problem  (BVP)  to  be  solved  at  each  time  step.  In  recent  years, 
the  MOLT  has  been  applied  as  a  numerical  method  to  solve  various  time-dependent  PDFs, 
with  initial  focus  on  parabolic  equations  [10,  19,  20,  6,  25].  The  work  [8]  extended  this 
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approach  to  the  second  order  wave  equation,  and  is  the  basis  of  the  present  work.  Further 
work  extended  this  approach  to  higher  dimensions  through  dimensional  splitting  [7]  and  to 
higher-order  through  the  Fourier  multiplier  approach  mentioned  in  the  previous  section  [9]. 
Based  on  the  work  in  [8,  7,  9],  we  give  here  an  overview  of  the  derivation  of  two  useful 
second-order  schemes  for  the  wave  equation,  their  extension  to  multiple  dimensions  through 
dimensional  splitting,  and  a  fast  numerical  algorithm  for  the  ID  problem.  In  Section  7,  we 
give  semi-discrete  Von  Neumann  analyses  and  dispersion  relations  for  these  semi-discrete 
schemes,  along  with  a  new  proof  that  the  fully  discrete  versions  are  unconditionally  stable, 
in  the  sense  of  Von  Neumann  analysis. 


3.2.1.  Dispersive  Semi-Discrete  Scheme.  As  in  the  previous  section,  we  substitute  the  fol¬ 
lowing  second-order  centered  discretization  into  the  wave  equation  ^Uu  —  Aw  =  f{x,t): 


u 


n 

tt 


-  2u^  + 
Af2 


Instead  of  proceding  by  the  Fourier  multiplier  approach  as  mentioned  previously,  we  apply 
an  averaging  technique  with  similar  results.  We  average  the  Laplacian  term  in  time 


An”  =  ^ A  +  2u^  -F  +  0(Af2) 

and  substitute  into  the  wave  equation.  Defining  a  =  2/{cAt)  and  manipulating  gives  the 
semi-discrete  scheme 


+  2m”  -1- 


Au^  +  ^f{x,D)  +  OiAt^). 


We  call  this  the  dispersive  semi-discrete  scheme,  since  all  terms  in  the  semi-discrete  dis¬ 
persion  relation  are  real- valued  (see  7.2). 


3.2.2.  Diffusive  Semi-Discrete  Scheme.  We  substitute  the  following  backward  difference  for¬ 
mula  (BDF)  discretization: 


u 


n+1 

tt 


2m”+i  -  5m”  +  4m”-^  -  m”-2 
Af2 


llAt^ 

12 


utttt{x,r]) 


into  the  wave  equation  -^Uu  —  Am  =  f{x,  t). 

Rearranging,  defining  a  =  y/2/{cAt)  and  dividing  by  gives  the  semi-discrete  scheme 


(-^A  +  l]  m”+^  =  1  (5m”  -  4m”-^  +  m”-2)  +  ^f(x,  t”+A  +  O(At^). 
\  J  2  ^ 


We  call  this  the  diffusive  semi-discrete  scheme,  due  to  the  presence  of  imaginary-valued 
terms  in  the  semi-discrete  dispersion  relation  (see  7.2).  The  (slight)  dissipation  of  this  method 
proved  useful  in  our  PIC  simulations,  and  we  use  this  scheme  (rather  than  the  dispersive 
scheme)  in  all  numerical  test  problems  in  this  work. 
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3.3.  Solution  of  the  Modified  Helmholtz  Equation.  As  seen  in  the  previous  section, 
both  the  dispersive  and  diffusive  semi-discrete  schemes  include  an  elliptic  BVP  to  be  solved 
at  each  time  step.  The  resulting  PDE  is  sometimes  called  the  modified  Helmholtz  equation 
[20].  In  contrast  to  the  usual  frequency-domain  Helmholtz  equation  Am  -|-  =  /,  the 

modified  Helmholz  equation  has  a  nonoscillatory  Green’s  function.  The  oscillation  in  the 
solution  of  the  wave  equation  is  supported  by  the  presence  of  multiple  time  levels  in  the  semi¬ 
discrete  equations.  Our  solution  strategy  is  to  use  the  well-known  technique  of  dimensional 
splitting  [29]  to  reduce  problems  in  multiple  dimensions  to  problems  in  one  dimension,  to 
which  we  apply  a  fast  integral  solution  method.  The  dimensionally-split  integral  solution 
naturally  leads  to  unconditionally  stable  numerical  schemes  with  computational  cost  and 
coding  complexity  comparable  to  explicit  schemes. 


3.3.1.  Dimensional  Splitting.  For  smooth  m(x,  ?/),  we  have 


~~  {dxx  +  dyy)  +  1  )  M  —  (  —  —  d^x  +  1  )  (  ~  ~  77^: 


a 


a 


4  '.^xxyy 


U 


So,  to  approximately  solve  (d^x  +  dyy)  +  l)  u  =  f  we  will  instead  solve 

- ^dxx  +  0  ( - ^dyy  -|-  M  =  /. 


where  we  have  introduced  the  splitting  error 


-KdxxyyU  =  0(c^At^ 


a 


To  solve 


——dxx  +  1  )  (  — 72^yy  ^  ^  ~  f 


a 


we  define  w  =  {—^dyy  +  l)  m  and  solve  the  following  one-dimensional  BVPs  “line  by  line”: 

w{x,  ■)  =  f{x,  •) 

+  K,y)  =  w{-,y) 

j 

where  appropriate  boundary  conditions  are  supplied  (see  [8,  9]).  To  facilitate  the  “line-by- 
line”  solution,  we  discretize  the  domain  with  a  uniform  Cartesian  grid. 


3.3.2.  Integral  Solution  Method  in  ID.  Consider  the  modified  Helmholtz  equation  in  one 
dimension,  -|-  1  j  w  =  f{x)  for  x  G  C  =  (a,  6),  with  appropriate  boundary  conditions 

imposed  at  x  =  a  and  x  =  b.  We  can  write  m(x)  =  u^{x)  +  u^{x),  where 

pb 

u^{x)  =  dx' 

u^{x)  = 


are  the  particular  and  homogeneous  solutions  and  where  A  and  B  depend  on  the  boundary 
conditions  imposed,  as  well  as  the  values  of  u^{a)  and  u^{b)  [8].  Our  fast  integral  solver 
consists  of  a  fast  convolution  algorithm  for  the  evaluation  of  the  particular  solution  u^, 
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along  with  appropriate  algorithms  for  evaluating  the  homogeneous  solution  ,  which  can 
be  viewed  as  boundary  correction  terms.  For  many  common  boundary  conditions,  the  co¬ 
efficients  A  and  B  for  the  boundary  correction  terms  can  be  found  by  applying  the  given 
boundary  conditions  and  solving  a  2  x  2  system  by  hand  for  A  and  B.  For  instance,  in  the  case 
of  homogeneous  Dirichlet  boundary  conditions  u{a)  =  u{b)  =  0  and  defining  7  = 
we  find 


A 

B 


7m^(6)  —  u^{a) 

1  —  72 

'yu^{a)  —  u^{h) 
1  —  72 


In  the  case  of  periodic  boundary  conditions,  we  find 


A 

B 


u^{h) 

1-7 

u^{a) 

1-7' 


With  some  further  consideration,  other  boundary  conditions  can  be  derived,  including 
outflow  (absorbing)  boundary  conditions  (for  the  underlying  wave  equation,  based  on  one¬ 
way  wave  equations).  For  further  details,  see  [8,  7,  9]. 


3.3.3.  Fast  Numerical  Evaluation  of  the  ID  Convolution  Operator.  The  convolution  operator 
giving  the  particular  solution  can  be  decomposed  as 

n  r°° 

f {r dx' 

px  poo 

=  2  /  dx'  +  -  dx' 

=:  /'-[/l(i) +  /«[/!  W 

for  a  given  function  /.  Meanwhile,  we  have  the  recursion  relations 

I^[f]{x)  =  e-^^^I^[f]{x  -  Ax) +  ^  f  /(a:')e-“(^-*')  dx' 

^  Jx—Ax 

n  T*-!-  /\  •T* 

/«[/] (a:)  =  e-“^"/^[/] {x  +  Ax)  +  ^  /  /(a:')e-“("'-")  dx'. 

Based  on  these  observations,  we  outline  the  fast  algorithm  for  the  numerical  evaluation 
of  this  convolution  operator  on  a  uniform  Cartesian  grid  developed  in  [8,  7,  9].  Consider 
the  convolution  operator  applied  to  a  function  /  supported  on  the  interval  (a,  6),  with  the 
convolution  also  to  be  evaluated  in  (a,  h).  The  interval  is  discretized  into  N  equal  subintervals 
of  length  Ax  =  {h  —  a)/N,  with  endpoints  xi  =  a,  =  Xj  +  Ax  for  j  =  1,  ...,N.  We  denote 
Ij  =  I[f]{xj),  =  I^[f]{xj)  and  if  =  /^[/](xj),  as  defined  above.  Further,  we  define  the 
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local  integrals 


2  Jxj.i 

pXj^l 

-  /  dx' 


j  =  2,...,iV  +  l 
j  =  l,...,N. 


Suppose  we  have  already  computed  these  Jj"  and  J^.  Setting  Jf  =  I§_^_i  =  0,  we  can  then 
perform  a  recursive  evaluation  of  the  convolution  integral  by  computing 

and  In+i-j  =  e~°‘^^lN+i-j+i  +  ^n+i-j  i  =  We  then  sum  Ij  =  for  each 

j  =  1,  ..N  +  1. 

The  method  for  evaluating  the  local  integrals  Jj"  and  depends  on  the  nature  of  the 
integrand.  For  the  particle  convolution  integral,  where  the  integrand  is  the  sum  of  particle 
shape  functions,  we  can  analytically  evaluate  the  local  integrals.  This  is  described  in  Section 
4.  For  general  integrands,  and  specifically  for  the  terms  involving  u  on  the  right  hand  sides 
of  the  semi-discrete  schemes,  we  numerically  evaluate  the  local  integrals  with  a  quadrature 
rule.  It  is  important  to  note  that  any  given  quadrature  rule  may  or  may  not  deliver  an 
accurate  and  stable  overall  scheme  for  the  wave  equation.  To  achieve  accuracy  and  stability, 
we  use  a  quadrature  rule  found  by  analytically  integrating  against  a  Lagrange  polynomial 
interpolant.  For  a  quadratic  interpolant,  leading  to  second-order  accurate  quadrature  in 
space,  we  obtain  the  following  approximations 


~  Pf{xj)  +  Qf{xj-i)  +  R{f{xj+i)  -  2f{xj)  -h  f{xj_i)) 
jf  ^  Pf{xj)  +  Qf{xj+i)  +  R{f{xj+i)  -  2f{xj)  +  f{xj-i)) 

where,  defining  u  =  aAx  and  d  =  e~'^ , 


P 

Q 

R 


1  - 


1  —  d 

V 


-d  + 


1-d 

V 


1  —  d  1  “I"  d 


2u 


Higher-order  spatial  accuracy  can  be  obtained  by  using  higher-order  accuracy  quadrature 
rules,  with  the  outcome  of  unconditional  stability  limiting  the  choice  of  quadrature  rules. 
Further  details  can  be  found  in  [8,  7,  9].  In  the  Appendix,  Von  Neumann  analyses  are  carried 
out  for  the  fully  discrete  diffuse  and  dispersive  schemes  in  one  and  multiple  dimensions,  and 
it  is  shown  that  they  are  both  unconditionally  stable. 


4.  Implementation  of  Particles 

4.1.  Particle  Weighting  Scheme.  In  our  PIC  methods,  the  charge  and  current  densities 
are  represented  as  the  sum  of  particle  shape  functions: 
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P(x,t)  = 


i=l 


X-Xp,i(t) 


J(x,t)  =  g»Vp,^>5' 


2=1 


Ax 

x-Xp^^(t) 

Ax 


where  iVp  is  the  number  of  particles,  Xp  j(t)  and  Vp_j(t)  are  the  position  and  velocity, 
respectively,  of  particle  z  at  time  t,  and  *S'(x)  is  a  particle  shape  function.  It  should  be 
emphasized  that  these  are  not  physical  particles,  but  rather  macro-  or  superparticles  that 
represent  a  discretization  of  the  PDF  [4],  We  do  not  directly  weight  particle  charge  and 
current  density  to  the  grid,  but  rather  we  analytically  evaluate  the  particle  convolution 
integrals  corresponding  to  these  source  terms  with  the  algorithm  described  below. 


4.2.  Method  for  Controlling  Divergence  Error.  It  is  well  known  that  electromagnetic 
PIC  methods  which  do  not  satisfy  a  discrete  form  of  Gauss’  law  through  their  held  solvers 
and  charge  and  current  weighting  schemes  will  suffer  severe  numerical  errors  related  to  charge 
conservation  [26] .  We  seek  to  develop  staggered  grid  approaches  to  computing  potentials  and 
helds  that  satisfy  discrete  analogues  of  Gauss’  law  V  ■  E  =  p/eo  and  the  identity  V  ■  B  =  0. 
The  divergence-free  condition  for  B  will  be  easily  satished  in  general,  as  the  magnetic  held 
will  be  calculated  as  a  hnite  diherence  curl  of  the  vector  potential.  In  order  to  numerically 
enforce  Gauss’  law,  we  seek  to  perform  an  elliptic  divergence  correction.  Future  work  will 
consider  the  alternative  of  hyperbolic  divergence  cleaning. 

We  now  give  the  mathematical  underpinning  of  our  elliptic  divergence  correction  technique. 
The  electric  held  is  calculated  as  E  =  —  V<F  —  Then  Gauss’  law  may  be  rewritten  as: 


p/eo  =  V  ■  E 


=  -A$  - 


a(V- A) 

m  ' 


Thus,  the  scalar  potential  satishes  the  Poisson  equation 


-At  =  p/e,  +  ^ 

Our  method  is  based  on  the  observation  that  if  this  Poisson  equation  is  suitably  discretized 
and  solved  on  a  staggered  grid  to  provide  the  scalar  potential  used  in  calculating  the  electric 
held,  then  the  electric  held  will  automatically  satisfy  a  discrete  form  of  Gauss’  law.  While  the 
exact  form  of  the  staggered  grid  will  depend  on  which  components  of  the  current  density  and 
the  electric  and  magnetic  helds  are  retained  in  a  given  model,  our  method  guarantees  exact 
discrete  divergence  relations  independently  of  the  charge  and  current  weighting  schemes  used, 
of  the  nature  of  the  solver  used  for  A,  and  of  the  gauge  condition  specihed. 
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4.3.  Particle  Equations  of  Motion.  In  our  PIC  methods,  the  approximation  of  the  evo¬ 
lution  of  the  Vlasov  equation  amounts  to  the  integration  of  the  equations  of  motion  of  the 
particles: 


dx 


p,i 


dt 

dvp^i 


dt 


where  a^  j(t)  is  the  acceleration  of  particle  i  at  time  t.  To  evolve  the  particle  equations 
of  motion,  we  use  standard  numerical  methods,  such  as  the  explicit  leapfrog  method,  and 
obtain  particle  accelerations  through  the  usual  interpolation  of  helds  from  grid  points  to 
particle  locations  [4].  Fields  are  calculated  as  hnite  differences  of  potentials  on  the  grid. 
As  these  aspects  are  standard  and  not  the  focus  of  the  present  work,  we  do  not  elaborate 
further. 


4.4.  Fast  Convolution  Algorithm  in  ID.  We  now  describe  the  algorithm  used  for  the 
fast  exact  evaluation  of  the  convolution  of  charge  and  current  density  source  terms  due  to 
particles.  It  has  two  main  steps.  There  is  a  local  deposit  step  and  then  a  recursive  sweep 
step.  This  basic  structure  is  the  same  in  all  dimensions.  For  definiteness,  we  describe  the 
application  of  the  algorithm  to  linear  particle  shapes  in  one  and  two  dimensions.  However,  it 
may  be  generalized  to  any  separable  particle  shapes  with  compact  support  in  any  dimension. 
Note  that  this  includes  many  widely  used  particle  shapes  in  PIC  algorithms,  namely  typical 
spline-based  particle  shapes  and  (suitably  cut-off)  Gaussian  particle  shapes.  For  the  case  of 
the  charge  density  integral,  the  particle  shape  function  Sp  below  is  replaced  by  qpSp  and  its 
contribution  summed  to  the  charge  density  integral,  and  for  the  case  of  the  current  density 
integral,  Sp  is  replaced  by  'VpQpSp  and  its  contribution  summed  to  the  current  density  integral. 


4.4.1.  Local  Deposit  Step  in  ID.  Consider  particle  p  located  in  the  cell  [xm,  Xm+\\  in  a  uniform 
grid  with  cell  length  Ax.  Let  Sp{x)  be  the  shape  function  of  the  particle.  Assume  that  the 
support  of  Sp  has  length  2rAx  for  some  integer  r.  The  local  deposit  step  then  consists  of 
analytically  evaluating  the  integrals 


for  j  =  m  —  r,  ...,m,  ...,m  -|-  r  and  for  each  particle  p,  and  summing  their  values  on  the 
grid. 

For  linear  particle  shapes  (corresponding  to  r  =  1),  we  have  Sp{x)  =  S{^^^^r),  where 


S{x) 


[l  —  X 

X 

F 

X 

Let  a  =  {xp  —  Xm)/Ax,  where  Xp  is  the  location  of  the  particle.  For  simplicity,  let  Xp  =  0 
and  Ax  =  1.  For  linear  particle  shapes,  we  then  have  the  situation  displayed  in  Figure  1. 
The  desired  integrals  are  then  easily  evaluated  for  linear  particle  shapes: 
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—a 


Jt^  =  a  j  (1  +  dx' 

=  (((l-a)a-l)  +  e"“e-“)/a 

Jm+i  =  a  [  (1  —  dx 

J  —a 

=  (e“"((a  —  l)a  +  1)  +  aa  +  1  —  2e“"e"“)/Q; 

=  a  [  (1  -  dx' 

Jl-a 

=  (e-"(-aa  -  1)  +  e-“e““)/a 

=  aj  (1  +  dx' 

=  (e““((a  —  1)q;  —  1)  +  e“"“)/Q; 

/I— a 

(1  -  dx' 

■a 

=  ((1  —  a)a  +  1  +  e“"(— ao;  +  1)  —  2e“"“)/Q; 

=  a  [  (1  -  dx' 

Jl-a 

=  (aa  —  1  +  e“"“)/a 


Note  that  just  one  evaluation  of  an  exponential  function  is  required  per  particle  (namely 
e"“).  To  account  for  arbitrary  Ax,  we  make  the  substitution  a  a  Ax  =  v. 

To  obtain  the  total  local  deposit,  we  simply  sum  the  particle  contributions  on  to  the  grid. 
Let  Np  be  the  total  number  of  particles.  The  algorithm  for  the  local  deposit  step  is  given 
by: 

Initialize  =  0  for  all  k. 

for  p  =  1  :  Ap  do 

Particle  p  located  in  cell  [xm,Xm+i\- 


Figure  1.  S{x)  =  1  -  |x|,  |x|  <  1,  S'(x)  =  0,  |x|  >  1 
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for  j  =  m  —  r:m  +  r  do 
Compute  jf+i, 

Deposit  +  jf/,,  Jf  =  Jf  +  jf  ^ 

end  for 
end  for 

Note  that  the  local  deposit  step  costs  0{Np)  operations. 

4.4.2.  Recursive  Sweep  Step  in  ID.  Once  we  have  performed  the  local  deposit  step,  we 
complete  the  evaluation  of  the  particle  integral  with  a  recursive  sweep  step.  Suppose  we 
have  N  gridpoints,  Xi,  ...,X]sf.  The  algorithm  for  the  recursive  sweep  step  is  given  by: 
Initialize  If'  =  I§  =  0 
for  j  =  1  :  —  1  do 

jL  _  jL  ip-ujL 

-'i+1  “  '^i+i  ^  ^  -'i 

tR  —  jR  _i_  jR 
J-N-j  —  '^N-j  w  e  -'AT-j+i 

end  for 

I  =  1^  +  1^ 

Note  that  the  recursive  sweep  step  costs  0{N)  operations. 


4.5.  Fast  Convolution  Algorithm  in  2D.  For  a  separable  particle  shape  S{x,y)  = 
Sx{x)Sy{y),  we  have 

4|4|S]|(x,j/)  =  /r|SJ(a:)  •  41SJ(j/) 

=  +  4“|SJ(:t))  '  (/f|S„](B)  +  /J'|S„](!/)) 

=  /i-|SJ(l)  .  /°|S„1(»)  +  ISJ(x)  .  4‘'|S,](»)  + 

+  /Jls.Kx)  ■  /°|S„1(!/)  +  4"|sj(x) .  4^[sj(9) 

This  is  suggestive  of  how  we  will  build  the  2D  algorithm. 


4.5.1.  Local  Deposit  Step  in  2D.  Consider  particle  p  centered  at  {xp,yp)  G 
[yn,yn+i],  with  separable  particle  shape  Sp{x,y)  =  where 


S{x) 


[l  —  X 

X 

F 

X 

Xr 


Xm+l]  X 


The  support  of  the  particle  shape  is  shown  in  Figure  2. 

In  the  local  deposit  step,  we  form  a  tensor  product  on  the  grid  as  suggested  by  the  above 
decomposition.  Note  that  a  total  of  12  local  integrals  must  be  evaluated  for  each  particle, 
then  summed  onto  the  grid  as  a  tensor  product. 

Initialize  =  0  for  all  j,  k. 

for  p  =  1  :  Ap  do 

Particle  p  located  in  cell  [xm,Xm+i\  x  [yn,yn+i\- 
for  j  =  m  —  r\  m  +  r,  k  =  n  —  r  ■.  n  +  r  do 
Compute 
Deposit: 

jLD  _  jLD  _i_  jL,p  jD,p 

"'j  +  l,fc+l  ~  "'j+l,fc+l  w  '^j+l  ■ 
tLU  _  jLU  I  jL,p  jU,p 
~  W  "'i+1  ■  "'fc 
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Figure  2.  The  support  of  a  linear  particle  shape  Sp{x^y)  in  2D. 


jRD  _  jRD  I  jRx  7 

'-',  ^4-1  “T  ■  Ji 


'j,k+l 


R,p  jD,p 
'k+1 


3 


j,k 

end  for 

end  for 


jRU  _  jRU  I  tR,P  tU,p 
'J-i  k  ~  ^  j^k  “T 


4.5.2.  Recursive  Sweep  Step  in  2D.  The  recursive  sweep  step  is  similar  to  the  ID  case,  and 
is  given  below, 
for  k  =  1  :  Ny  do 

Initialize  =  0. 

for  j  =  1  :  iVa;  —  1  do 


tLD 


=  XD  p  +  e-"/, 


—u  tLD 


jLU  ^  JLU  ,  „ 
^j+l,k  '^j+l,k  '  ^ 


j,k 


tRD  _  jRD  I  „—u  tRD 

^N,-j,k  -  ^N,-j,k  +  e  J7v,-j+l,fc 

tRU  _  jRU  I  „—u  tRU 

^N^-j,k  —  ^N^-j,k  +  ^ 

end  for 
end  for 

jD  _  jLD  _j_  jRD 
jU  ^  jLU  _|_  pRU 

for  j  =  1  :  iVa;  do 

Initialize  iR  =  =  0. 

ja  JD^y 

for  A;  =  1  :  iV,,  —  1  do 


y 

=  jp 


j,k+l  '  ^  -^j,k 


end  for 
end  for 

I  =  lD  +  iu 


jU  _  jU 

^j,Ny-k  —  ■'^j,Ny-k  ^  ^  ^j,Ny-k  +  l 
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Both  ID  and  2D  overall  algorithms  cost  0{Np  +  N)  operations,  where  N  is  the  total 
nnmber  of  gridpoints.  Since  in  a  typical  PIC  simnlation,  Np  »  N,  the  cost  of  the  overall 
algorithm  is  dominated  by  the  local  deposit  step. 


4.6.  Particle  Boundary  Conditions.  In  dealing  with  bonndaries,  two  types  of  consider¬ 
ations  mnst  be  made.  First,  we  mnst  determine  what  to  do  in  the  integration  of  bonndary 
particles,  for  which  the  snpport  of  the  shape  fnnction  extends  ontside  of  the  domain.  This 
will  be  dependent  npon  the  type  of  bonndary  condition.  For  periodic  bonndary  conditions, 
we  can  simply  extend  the  particle  shape  fnnction  periodically  and  proceed  to  integrate.  For 
Dirichlet  bonndary  conditions,  we  extend  the  integration  domain  to  inclnde  ghost  points  jnst 
beyond  the  bonndary  to  which  the  bonndary  particles  are  weighted. 

Second,  we  mnst  ensnre  that  the  particle  convolntion  integral  is  consistent  with  the  bonnd¬ 
ary  conditions  on  the  wave  fnnction.  This  is  easily  handled  throngh  the  nsnal  bonndary 
correction  terms  in  one-dimension,  and  can  be  extended  to  the  dimensionally-split  mnltidi- 
mensional  case. 


5.  Numerical  Results 

5.1.  Electrostatic  Test  Problems.  We  first  consider  three  standard  periodic  electrostatic 
test  problems  in  ID  and  2D,  then  a  ID  bonnded  plasma  problem,  the  simnlation  of  sheath 
formation.  In  the  first  three  test  problems,  electrons  are  loaded  from  a  pertnrbed  initial 
distribntion  of  the  form 


fe{x,v,t  =  0)  =  fe{v)  (^1  +  e sin  (5) 

where  is  the  length  of  the  domain,  e  is  the  amplitnde  of  pertnrbation,  and  fe{v)  is  the 
initial  velocity  distribntion.  In  the  2D  case,  simnlations  are  taken  to  be  nniform  in  the  y- 
direction.  We  normalize  qnantities  according  to  the  nondimensionalization  presented  in  the 
Section  7.1.1.  In  particnlar,  we  normalize  time  qnantities  to  the  inverse  plasma  freqnency 
ojp^.  We  will  consider  a  periodic  domain  with  a  nniform  nentralizing  backgronnd  charge, 
and  fnrther  we  set  the  speed  of  light  c  =  100.  In  these  problems,  we  see  good  performance 
even  at  large  CFLs,  since  the  physics  is  dominated  by  the  low  freqnency  spatial  modes. 

5.1.1.  Cold  Plasma  Langmuir  Wave.  We  consider  a  cold  plasma  Langmnir  wave  [4]  with 
/e(n)  =  d{v).  Electrons  are  pertnrbed  away  from  a  nniformly  distribnted,  motionless  state 
against  a  static,  nniform  nentralizing  backgronnd  charge  distribntion.  The  resnlting  separa¬ 
tion  of  charge  prodnces  cold  plasma  oscillation.  In  the  ID  simnlation,  we  set  =  2n  and 
e  =  0.1.  We  nse  a  100  cell  grid  and  take  At  =  0.1,  and  we  nse  Np  =  3600  particles.  In 
the  2D  simnlation,  we  set  =  Ly  =  2tt  and  the  pertnrbation  strength  e  =  0.1.  We  nse  a 
100  X  100  grid  and  again  take  At  =  0.1,  and  we  nse  Np  =  360000  particles.  In  both  ID  and 
2D  cases,  the  CFL  nnmber  nsed  is  cAt/Ax  ~  159,  mnch  larger  than  what  wonld  be  allowed 
by  an  explicit  method.  The  oscillation  in  the  potential  energy  is  plotted  and  compared  to 
the  prediction  of  linear  theory  in  Fignre  3;  we  see  that  the  plasma  freqnency  is  accnrately 
reprodnced. 
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Figure  3.  Potential  energy  in  cold  plasma  oscillation.  Green  is  the  ID  nu¬ 
merical  result,  blue  is  the  2D  numerical  result,  and  red  is  the  prediction  of 
linear  theory.  We  see  the  plasma  frequency  is  accurately  reproduced  in  our 
simulations. 


5.1.2.  Two  Stream  Instability.  We  consider  the  two  stream  instability  with  fe{v)  =  6{v  — 
Weam)  +  5(n  -f  fbeam)-  Two  counterstreamiug  beams  of  electrons  are  perturbed  away  from 
a  uniformly  distributed  state  against  a  static,  uniform  neutralizing  background  charge  dis¬ 
tribution.  The  beams  interact  and  “roll  up”  in  phase  space,  causing  some  of  the  particles’ 
kinetic  energy  to  be  transformed  into  potential  energy  stored  in  the  electric  held.  According 
to  the  dispersion  relation  for  the  two  stream  instability  from  linear  theory  [4],  we  have 


ui 


2u^(ujI 


beam 


2ujt)  =  0 


(6) 


which  gives  the  greatest  growth  rate,  7  ~  0.3535,  for  k  ~  3.06.  We  therefore  scale  the 
domain  to  this  value  of  A;,  and  take  =  27r/3.06.  We  take  the  beam  velocity  Ubeam  =  1 
and  the  perturbation  strength  e  =  0.0005.  In  our  ID  simulation,  we  use  a  100  cell  grid 
with  At  =  0.1,  and  we  use  Np  =  1000  particles.  In  our  2D  simulation,  we  use  a  100  x  100 
grid  with  At  =  0.1,  and  we  use  Np  =  1000000  particles.  This  results  in  a  CFL  number  of 
cAt/Ax  ~  68.  We  run  the  simulations  for  1000  time  steps.  The  growth  of  the  k  =  3.06 
mode  of  the  electric  held  is  shown  in  Figure  4  for  the  ID  and  2D  cases,  and  agrees  with 
the  rate  from  linear  theory.  In  the  nonlinear  saturation  stage,  we  see  a  slight  discrepancy 
between  the  ID  and  2D  results,  probably  due  to  the  accumulation  of  numerical  error.  We 
also  show  selected  phase  space  plots  in  Figure  5,  where  we  see  the  expected  “rolling  up”  of 
the  two  beams.  Resolution  is  limited  by  the  number  of  particles. 


5.1.3.  Landau  Damping.  We  consider  Landau  damping  of  Langmuir  waves  in  a  warm  plasma, 
with  fe{v)  taken  to  be  Maxwellian.  Warm  electrons,  following  a  Maxwellian  velocity  distri¬ 
bution,  are  perturbed  away  from  a  unifom  distribution  against  a  static,  uniform  neutralizing 
background  charge  distribution.  Potential  energy  from  the  electric  held  is  transformed  into 
kinetic  energy  of  particles.  The  dispersion  relation  from  linear  theory  in  this  case  gives  a 
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Figure  4.  Growth  of  the  mode  with  maximum  growth  rate  in  the  two  stream 
instability,  corresponding  to  k  =  3.06.  Green  is  the  ID  numerical  result,  blue 
is  the  2D  numerical  result  (measured  along  the  central  y  =  0  slice),  and  red  is 
the  prediction  of  linear  theory.  We  see  the  correct  growth  rate  is  reproduced 
in  our  simulations. 


decay  rate  of  7  ~  0.154  for  the  k  =  0.5  mode  [4].  We  take  =  47r,  electron  thermal 
velocity  Utherm  =  1  and  perturbation  strength  e  =  0.1.  In  our  ID  simulation,  we  use  a  100 
cell  grid  and  take  At  =  0.1,  and  we  use  Np  =  1000000  particles.  In  our  2D  simulation,  we 
use  a  100  x  100  grid  and  take  At  =  0.1,  and  we  use  Np  =  9000000  particles.  We  run  the 
simulations  for  300  time  steps.  The  decay  of  the  k  =  0.5  mode  of  the  electric  held  in  the  ID 
and  2D  simulations  is  shown  in  Figure  6,  and  agrees  with  the  rate  from  linear  theory.  As  in 
the  two  stream  instability  example,  there  is  a  discrepancy  between  the  ID  and  2D  results  at 
later  times,  again  likely  due  to  the  accumulation  of  numerical  errors. 

5.1.4.  Sheath  Formation  in  a  Bounded  ID  Plasma.  We  present  the  simulation  of  sheath 
formation  in  a  bounded  ID  plasma,  following  the  model  described  in  [32].  In  contrast  to  the 
previous  problems,  this  simulation  incorporates  both  mobile  electrons  and  ions.  Electrons 
and  ions  are  initialized  from  Maxwellian  distributions  and  uniformly  spatially  distributed  in 
a  bounded  domain.  The  left  boundary  is  a  symmetry  plane,  and  so  we  impose  Neumann 
boundary  conditions  on  the  potential,  and  rehux  boundary  conditions  on  particles,  as  in 
[32].  The  right  boundary  is  a  conductor  that  collects  charged  particles.  When  particles  hit 
the  right  boundary,  they  are  removed  from  the  simulation.  Since  electrons  have  a  higher 
average  velocity  than  ions,  they  have  a  greater  flux  on  the  collector  and  become  depleted 
near  the  right  boundary,  where  the  difference  between  ion  and  electron  densities  leads  to  a 
collector  sheath  region,  where  the  potential  changes  from  the  interior  value  to  the  wall  value, 
which  has  the  effect  of  repelling  electrons  away  from  the  right  wall.  Electrons  and  ions  are 
replenished  by  a  particle  source  region  near  the  left  boundary,  where  electrons  and  ions  are 
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injected  uniformly  in  the  region  from  a  Maxwellian  distribution  at  a  fixed  rate  per  time  step. 
We  take  mj/rue  =  100,  Tsrc,i/Tsrc,e  =  1)  and  we  set  Vth,e  =  1,  and  set  =  20  (in  Debye 
lengths).  We  use  a  100  cell  grid  and  take  At  =  0.1,  which  gives  a  CFL  number  of  50.  We  run 
our  simulation  for  8000  time  steps,  up  to  3.6  thermal-ion  transit  times.  In  Figure  7  we  see 
the  result  of  the  simulation.  In  Figure  7a,  we  see  the  profile  of  the  potential,  which  has  the 
right  qualitative  features,  including  a  collector  sheath  region  that  is  several  Debye  lengths 
wide.  In  Figure  7b,  we  see  the  net  electron  and  ion  counts,  along  with  the  injection  rate. 
The  difference  between  the  electron  and  ion  counts  reflects  the  difference  between  electron 
and  ion  densities  in  the  collector  sheath  region. 


5.2.  Electromagnetic  Test  Problems. 


5.2.1.  Bennett  Pinch  Problem.  We  present  the  application  of  our  PIC  method  to  the  Bennett 
pinch  [3],  an  effect  related  to  the  magnetic  confinement  of  a  beam  of  charged  particles.  A 
beam  of  charged  particles  induces  a  solenoidal  magnetic  field  around  the  beam.  Particles 
near  the  edge  of  the  beam  move  orthogonally  to  these  field  lines  at  the  beam  drift  velocity, 
causing  the  particles  to  be  accelerated  towards  the  center  of  the  beam,  in  effect  confining 
particles  in  the  beam.  An  appropriate  choice  of  parameters  leads  to  a  stationary  steady  state, 
uniform  along  the  axis  of  the  beam.  A  well-known  magnetohydrodynamic  (MHD)  model  of 
a  stationary  steady  state  gives  explicit  formulas  for  the  beam  density  and  the  magnetic  field 
[5],  and  provides  a  basis  for  the  validation  of  our  numerical  method.  Moreover,  it  is  a  first 
step  toward  applying  our  method  to  more  physically  interesting  beam  instability  problems 
in  three  dimensions. 

Our  PIC  simulation  of  the  Bennett  pinch  is  two-dimensional  in  physical  space,  and  three- 
dimensional  in  velocity  space.  The  particle  beam  is  considered  uniform  along  the  axis  of  the 
beam,  which  we  take  to  be  the  2;-direction,  which  reduces  the  physical  dimensions  to  two. 
Electrons  drift  in  the  2;-direction  with  a  uniform  average  beam  drift  velocity  Uf,,  and  this 
motion  induces  a  confining  magnetic  field  with  only  x-  and  ^/-components.  A  stationary  ion 
background  distribution  enforces  quasineutrality  in  the  beam,  with  any  separation  of  charge 
producing  an  electric  field  with  only  x-  and  //-components  acting  as  a  restoring  force.  The 
electrons  are  assumed  to  follow  a  Maxwellian  distribution  with  thermal  velocity  Vth-  We  take 
Vb/vth  =  100  and  c/vth  =  1000,  where  c  is  the  speed  of  light.  The  ions  are  considered  cold 

m  =  0). 

Since  the  beam  drift  velocity  is  taken  to  be  much  larger  than  the  (transverse)  thermal 
velocity,  and  further,  the  transverse  velocities  follow  a  Maxwellian  distribution  and  so  should 
not  generate  any  net  currents,  we  neglect  the  x-  and  //-components  of  the  current  density 
(and  so  also  of  A).  In  the  true  solution,  ^  =  0  and  =  0,  so  we  neglect 
Hence,  we  actually  only  solve  two  wave  equation,  one  for  A^,  obtaining  only  transverse 
magnetic  field  components,  and  By  =  —  and  one  for  d),  obtaining  only 

transverse  electric  field  components  and  Ey  =  —  Thus,  the  Poisson  equation 

satisfied  by  the  scalar  potential  is  —  A<F  =  p/cq.  We  discretize  our  domain  with  a  staggered 
grid,  one  cell  of  which  is  shown  in  Figure  8. 

The  scalar  potential  is  calculated  from  the  standard  5-point  finite  difference  Laplacian, 
and  satisfies  the  equation 
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The  electric  and  magnetic  fields  are  calculated  on  the  staggered  grid  by  finite  differences 


/7'*4+1/2  _  _ 
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Ax 
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4*J+l  _  4*J 

Dij+l/2  _ 

"  Ax 

,  _  /iti 

D*+l/2,i  ^ 

^  Ax 

The  electric  field  then  satisifies  the  following  discrete  analogue  of  Gauss’  law; 


[V  ■  = 


pi+l/2,j  jpi,j+l/2  Tpij  — 1/2 


'■  ’  ■'  Ax  '  Ax 

1  //  /  <|)*j  _ 

Ax  \\  Ax  /  \  Ax  / 

/  <|)* J+l  _  <|)*4  \  /  <|)ij  _  \ 

\  Ax  )  \  Ax  ) ) 

Ax^ 

=  p'V^o. 

The  magnetic  field  satisfies  the  following  discrete  analogue  of  the  divergence  free  condition; 


ji+lj+l/2  j^i,j+l/2  T^i+1/2J+1  j^i+l/2,j 


[V  ■  =  — _ ~  _ 

Ax 

Ax  \\  Ax 

/  ^i+ij+i  _  ^*4+1  \ 

^  (  A^^^  )  “ 


-  Ai’A 

A^^^  )  ^ 


_  ^i  j 
z  z 

Ax 


All  computational  boundaries  in  this  problem  are  outflow  boundaries.  In  order  to  supply 
the  finite  difference  Poisson  solver  with  suitable  boundary  values,  the  wave  solver  is  applied 
with  outflow  boundaries  conditions  to  evolve  the  wave  potential  alongside  the  Poisson 
potential  <h.  The  boundary  values  from  <hvv  are  then  supplied  to  the  Poisson  solver  to  use 
in  calculating  $.  Once  the  wave  solver  reaches  steady  state,  <hvv  arid  $  differ  only  by  0.1% 
relative  error,  however,  the  wave  potential  gives  a  discrete  electric  field  with  divergence  error 
on  the  order  of  10“^,  while  the  Poisson  potential  gives  a  discrete  electric  field  with  divergence 
error  on  the  order  of  machine  epsilon  10“^®. 
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Like  in  the  other  test  problems,  we  use  the  diffusive  version  of  wave  solver.  Particle 
velocities  in  all  three  directions  are  updated  with  the  nonrelativistic  Boris  push  [4],  Particles 
are  initialized  according  to  the  MHD  steady  state  (according  to  the  theoretical  spatial  density 
profile  and  the  corresponding  Maxwellian  distribution  in  velocity  space)  and  held  fixed  while 
the  field  solver  is  stepped  to  an  approximate  steady  state,  after  which  the  particle  push  is 
turned  on.  The  simulation  is  run  to  a  final  time  of  Rb/{2vth)  (plus  startup  time),  where 
Rb  is  an  effective  beam  radius  and  Vth  is  the  thermal  velocity,  at  which  time  there  would 
be  substantial  spreading  of  the  beam  in  absence  of  the  confinement  effect.  We  choose  Rb 
such  that  99%  of  the  particles  in  the  theoretical  beam  are  within  this  radius.  In  loading 
particles,  the  beam  is  cut  off  at  radius  Rb  (no  particles  are  loaded  outside  of  this  radius). 
The  computational  domain  is  taken  to  be  a  ARb  x  4Rb  square  centered  on  the  beam  axis.  The 
computational  domain  is  truncated  with  outflow  boundary  conditions.  Particles  exiting  the 
boundary  of  the  computational  domain  are  reinjected  into  the  beam  to  maintain  constant 
total  current.  However,  since  most  particles  should  be  confined  within  the  beam,  such 
boundary  crossings  should  be  rare. 

Numerical  results  for  the  Bennett  pinch  are  given  in  Figure  9.  In  order  to  resolve  large 
gradients  near  the  center  of  the  beam,  we  use  a  500-cell  by  500-cell  grid  and  a  CFL  num¬ 
ber  of  cAt/Ax  =  3  (except  in  Figure  9c  as  noted)  and  500,000  electron  particles.  Fi¬ 
nal  numerical  solutions  are  shown  at  the  final  time,  after  334  start  up  time  steps  and 
20,834  PIC  time  steps.  (The  final  time  is  approx.  22  plasma  periods,  and  the  diameter 
of  the  beam  is  approx.  280  Debye  lengths.)  In  Figure  9a,  we  see  good  agreement  be¬ 
tween  the  numerical  electron  density  and  MHD  theory.  The  inset  zoomed  portion  shows 
a  slight  discrepency  at  the  peak  of  the  beam,  due  to  statistical  fluctuation  caused  by  the 
hnite  number  of  particles.  In  Figure  9b,  we  see  the  time  histories  of  the  potential  en¬ 
ergy,  calculated  as  ^(H^  ^  -|-  H^^^)  -|-  -\-  E^^^)  where  the  sum  is  over  grid 


points  j  (with  er  and  /ir  defined  as  in  Section  7.1.2),  and  the  kinetic  energy,  calculated  as 
+  '^y,i  +  —  VbY)  where  the  sum  is  over  electron  particles  i.  We  see  good 

energy  conservation,  despite  the  slight  dissipation  of  the  diffusive  scheme.  The  initial  spike 
in  the  potential  energy  is  the  result  of  transient  waves,  arising  due  to  the  beam  turning  on, 
and  flowing  out  of  the  domain  as  the  solution  is  stepped  to  a  steady  state.  In  Figure  9c, 
we  see  the  result  of  refinement  in  At,  keeping  Ax  fixed,  showing  a  profile  of  the  azimuthal 
magnetic  field  Bg  along  the  central  y  =  slice  for  CFL  numbers  of  3,  10  and  20,  along  with 
MHD  theory.  We  observe  approximate  second-order  convergence  in  At,  as  expected  (a  more 
robust  convergence  study  is  confounded  by  the  slow  convergence  in  particle  number  in  PIC 
methods).  Outside  of  the  beam  radius  Rb  =  1,  there  is  error  associated  with  the  finite  cut-off 
radius  of  the  beam  (the  theoretical  beam  density  decays  only  algebraically).  In  Figure  9d, 
we  see  the  numerical  error  of  the  azimuthal  magnetic  field  Bq  (with  a  CFL  number  of  of 
3)  normalized  by  the  peak  value  of  the  magnetic  field,  and  we  see  that  there  is  a  geometric 
pattern  to  the  numerical  error,  characteristic  of  the  dimensionally-split  method.  In  addition 
to  this  splitting  error,  the  total  error  is  contributed  to  by  errors  associated  with  the  spatial 
quadrature  and  the  finite  differences  used  to  calculate  the  magnetic  field  (likely  contributing 
to  the  large  error  at  the  center  of  the  beam  due  to  large  gradients  there)  and  with  the  finite 
beam  cut-off  radius  and  the  outflow  boundary  condition  (contributing  most  strongly  near 
the  boundary  of  the  computational  domain).  These  results  show  that  our  method  can  indeed 
simulate  a  basic  electromagnetic  plasma  phenomenon  with  a  CFL  number  larger  than  what 
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is  allowed  by  typical  explicit  schemes.  The  CFL  number  used  in  this  problem  is  limited  by 
the  accuracy  of  the  second-order  wave  solver.  A  higher-order  wave  solver,  such  as  those  in  [9] , 
would  allow  for  a  larger  usable  time  step  size,  and  will  be  the  subject  of  further  investigation. 


5.2.2.  Mardahl  Beam  Problem.  We  apply  our  method  to  the  beam  problem  proposed  in  [26] 
as  a  diagnostic  for  the  effects  of  divergence  error. 

In  the  Mardahl  beam  problem,  we  have  currents  in  the  plane  of  simulation  only,  and  so 
we  retain  the  x-  and  y-components  of  the  vector  potential,  and  Ay,  along  with  the  scalar 
potential  <h  in  the  model.  We  retain  the  electric  held  components  ~  and 

Ey  =  and  the  magnetic  held  component  The  Poisson  equation 

satished  by  the  scalar  potential  is 


-A$  =  p/eo  + 


m 


dA^  dAy\ 

dx  dy  ) 


We  discretize  our  domain  with  a  staggered  grid,  one  cell  of  which  is  shown  in  Figure  10. 
Denoting  by  D^t  a  (linear)  hnite  diherence  discretization  of  the  time  derivative  operator 
d/dt,  the  scalar  potential  satishes  the  equation 


Ax'^ 


=  P^Veo+ 

4*+l/2,j 

_  ,  Jl.'T' 

Dm 


A 


*-1/2, j 


Ax 


+ 


a: 


iJ+l/2 

y 


a: 


y 


Ax 


The  electric  and  magnetic  helds  are  calculated  on  the  staggered  grid  by  hnite  diherences 
as 


rpi+l/2,j 

■‘-^X 

rpi,j+ll2 

y 

gi+l/2,j+l/2 


4)*+i,i  _ 

Ax 

<|)*,i+i  _  4)*,i 
Ay 

K"''  - 

Ax 


DA,(Ay*^P 


Ay 


The  electric  held  then  satisihes  the  following  discrete  analogue  of  Gauss’  law: 
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[V  ■  E]*’^'  = 


e: 


|*+l/2,i  _  77t*J+1/2  jpi,j—lf2 


1 

Ax 


.  - h 

Ax 


y 


-E, 


Ax 
Ax 


Ay 


<i)*j  _  <|)*-ij 


Ax 

<|)*J  _ 


-  I  + 


Ax 


<j)*+ij  _|_  <j)*j+i  _j_  _|_  _  4$*’-? 

Ax'^ 


D 


At 


*+l/2j  _  Ai-l/2,j  Ai,j+l/2  _  Ai,j-l/2' 

X  -^X  , 

- h  — 


Ax  Ax 

=  p"V^o. 

where  we  have  used  the  linearity  of  D/^f 


6.  Conclusions  and  Future  Work 

In  this  work,  we  have  described  a  PIC  method  that  uses  an  unconditionally  stable  wave 
equation  solver  to  eliminate  the  CFL  restriction  on  the  ratio  of  the  time  step  size  to  the 
spatial  step  size,  typical  of  explicit  methods,  while  retaining  computational  cost  and  code 
complexity  comparable  to  such  explicit  methods.  Our  numerical  results  show  that  we  can 
apply  our  method  to  problems  of  plasma  physics  using  a  time  step  size  larger  than  what 
would  be  allowed  by  a  typical  explicit  field  solver.  We  have  seen  that  the  usable  time  step 
size  can  be  limited  by  the  numerical  accuracy  of  the  method  when  there  are  large  gradients 
(high-frequency  content)  in  the  solution.  We  implement  a  staggered  grid  approach  to  give  an 
electromagnetic  PIC  method  that  satisfies  exactly  a  discrete  form  of  Gauss’  law.  Future  work 
will  center  on  making  use  of  the  implicit  wave  solvers  ability  to  handle  complex  boundary 
geometries  without  the  use  of  a  staircasing  approximation.  We  will  investigate  the  use  of 
higher-order  methods,  such  as  given  in  [9],  in  our  PIC  method  in  order  to  increase  the 
maximum  usable  time  step  size  (we  have  already  found  a  4th  order  Newmark  scheme  with 
slight  dissipation  that  may  prove  useful).  A  further  course  of  action  will  be  to  implement  a 
boundary  integral  treecode  (BIT)  solution  to  solve  the  modified  Helmholtz  equations  in  the 
semi-discrete  schemes,  such  as  in  [24],  rather  than  use  dimensional  splitting. 

7.  Appendix 

7.1.  Nondimensionalization  and  Asymptotic  Analysis.  We  provide  the  normaliza¬ 
tions  used  in  the  test  problems  presented  in  this  work  in  both  electrostatic  and  electromag¬ 
netic  cases.  In  the  electrostatic  case,  we  argue  by  formal  asymptotic  analysis  and  classical 
solution  formulas  that  the  Vlasov- Wave  system  agrees  with  the  Vlasov-Poisson  system  in 
the  relevant  scaling  limit. 

7.1.1.  Electrostatic  Case.  Here  we  give  the  normalization  used  for  the  electrostatic  test  prob¬ 
lems  in  Section  5.1.  Consider  the  following  change  of  variables: 

f  =  Ef,  i=Tt,  X  =  Lx,  V  =  Vv,  p  =  qNp,  l>  = 

22 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


applied  to  the  system 


^  +  v-Vj-^Vl>-Vj=0 
at  m 


1  (92$ 


A$  =  p/£o,  p{x,t)  =  q  dir 


Assuming  the  scalings, 


y  T 

T’ 


meo  qNL"^  {eomY^'^ 

2  ^  '^p  ’  *i’o  =  — '  y  = 


Nq 


^0 


q^L° 


which  are  the  natural  scalings  in  the  electrostatic  limit,  we  obtain; 


|^+v-Vx/-V$-Vv/  =  0 


0.92$ 


-  A$ 


p,  p(x,t) 


/(x,v,t)  dv 


where  e  =  ^  ^  (not  to  be  confused  with  the  electric  permittivity  Sq)-  Note  that  1/e  is 

the  speed  of  propagation  of  waves  in  the  potential  in  this  normalization,  and  becomes  large 
when  e  is  small,  that  is,  when  the  characteristic  particle  velocities  are  small  compared  to  the 
speed  of  light. 

Assume  the  following  formal  asymptotic  expansions; 


/  —  /o  +  e/i  +  e^/2  + 
p  =  [  fdv 


=  /  fod'v  +  €  /idv  +  eM  /2dv  H - 

Jv  J  w  J  w 

=  Po  +  epi  +  e2p2  +  ■■■ 

$  =  $0  +  e$i  +  e2$2  H - 
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Collecting  in  orders  of  e: 


0(1): 


0(e): 


O(e^),  k>2- 


^  +  v- Vx/o-V$o- Vv/o  =  0 

-A<l'o  =  po,  Po=  fod\ 

J  V 

^  +  V  ■  Vx/i  -  V<l>0  ■  Vv/i  -  V$1  ■  Vv/o  =  0 
at 

-A$i  =  pi,  Pi  =  /  /idv 


dfk 

dt 


+  V  ■  Vx/fc  -  ^  V<l>j  ■  Vv/fc-j  =  0 

j=o 


Acf>  f  f  W 

-A<I)fc  =  Pfc - pk=  Ik  ctv 


We  note  that  the  leading  order  is  precisely  the  Vlasov-Poisson  system  (nondimensionalized 
under  the  same  scalings).  This  formal  computation  suggests  that  our  model  will  agree  with 
the  electrostatic  model  to  0(e  =  V/c)  when  particle  velocities  are  small  compared  to  the 
speed  of  light.  This  can  be  considered  as  a  consequence  of  the  strong  Huygens’  principle  in 
odd  spatial  dimensions,  and  of  a  weaker  decay  property  that  holds  in  even  spatial  dimensions, 
which  can  be  deduced  from  classical  solution  formulas  [15,  30].  Consider  the  Cauchy  problem. 


1  d‘^u 


dt^ 


Au  =  f{x,t) 


{x,  t)  G  X  (0,  oo) 


u{x,  0)  =  0 
ut{x,0)  =  0 


X  G 
X  G 


for  d  =  2,3  and  /  sufficiently  smooth  with  compact  support,  for  which  classical  explicit 
solution  formulas  exist.  We  consider  the  case  of  /  having  compact  support  in  B{0,R)  x 

d 

(0,T)  C  M'’*  X  (0,  cxd),  where  =  {x  G  M'’*|  X)  Classical  solution  formulas 

i=i 

imply  that  for  x  G  5(0,  R)  and  t  >  T  +  2R/cT,  we  have 

■u(x,  t)  =  0{{ct)~^) 
u{x,  t)  =  0 


d  =  2 
d  =  3 


As  a  generalization  of  this,  for  any  sufficiently  smooth  f{-,t)  supported  in  B{0,R)  for  all 
f  >  0  with  =  /r(-)  for  all  f  >  T,  it  is  again  easily  argued  that  for  x  G  B{0,R)  and 

t  >  T  +  2R/cT,  we  have 


u{x,t)  =  u\,{x)  +  0{{ct)  ^) 
u{x,  t)  =  Mp(x) 


d  =  2 
d  =  3 


where  Mp(x)  is  the  classical  integral  solution  of  the  Poisson  equation  —Au%  =  fr  in 


dimension  d. 
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The  convergence  of  solutions  to  the  Vlasov-Maxwell  system  to  those  of  the  Vlasov-Poisson 
system  has  been  rigorously  considered  in  works  such  as  [2,  34],  It  may  be  possible  to  apply 
similar  techniques  to  rigorously  study  the  convergence  of  solutions  of  our  Vlasov-Wave  model 
to  those  of  the  Vlasov-Poisson  system,  but  this  is  outside  of  the  scope  of  the  present  work. 


7.1.2.  Electromagnetic  Case.  Here  we  give  the  normalization  used  for  the  electromagnetic 
test  problem  in  Section  5.2.1.  Consider  the  following  change  of  variables: 

f  =  Ff,  i=Tt,  x  =  Lx,  V  =  Vv, 

p  =  qNp,  J^  =  qVNJ^, 

$  =  $0$,  A  =  AqA, 

applied  to  the  system 

+  V  ■  Vx/  +  -  (- Vxl>  -  V  X  (Vx  X  (0, 0,  A,))]  ■  Vv/  =  0 

1  92$ 


c2  dC 

1 


c2  dC 
Assuming  the  scalings, 


-  A<h  =  p/so, 
AA^  PqJzi 


p(x,t)  =q  /(x,v,t)civ 
J  V 

J^{5t,i)  =  q  /  h^/(x,  v,t)dv 


T/  _  ^  m  _  -1 

^  T’  V  Ag2  ’ 

^  qNL^en  ,  poVqNL^  /  2 

=  - ,  An  =  - ,  ErPr  =  V  /c  , 


F  = 


^0 

{eoiriYF 

qdLd 


where  Er  and  pr  are  dimensionless  parameters  introduced  to  enforce  the  Lorenz  gauge 
condition,  we  obtain: 


+  V  ■  Vx/  +  £R(-V<h  +  V  X  (V  X  (0, 0,  A,))  ■  Vv/  =  0 


where  e  =  ^  =  7, 
5.2.1,  we  choose  Er 


^  -  Ad>  =  p/er,  p(x,  t)=  /(x,  V,  t)  dv 


'  dC 


AA^  =  PrJ^,  J^(x,t)  =  /  n^/(x,  v,t) 

J  V 


dv 


as  in  the  electrostatic  case.  For  the  Bennett  pinch  problem  in  Section 

=  1. 


7.2.  Semi-Discrete  Von  Neumann  Analysis.  In  this  section,  we  provide  semi-discrete 
Von  Neumann  stability  analyses  and  dispersion  relations  for  the  semi-discrete  schemes  de¬ 
rived  in  3.  These  build  on  similar  analyses  for  related,  but  different,  schemes  given  in  [8]. 
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7.2.1.  Dijfusive  Scheme.  Substituting  the  ansatz 

^  (5u^  -  4u^-^  +  , 

we  obtain  a  polynomial 


-  4A^  +  5A  -  2(1  +  z^)  =0 

where  A  =  =  |A;|/q;.  The  three  roots  tell  us  about  possible  modes 


A  necessary  condition  for  stability  is  A  >  1  for  all  roots. 

The  hrst  root, 

Ai  =  —  +  —  ^27^:^  +  3y/3y/27z^  +  2z‘^ 

1 

3  {27 z^  +  3V3V27z^T^  +  l)  ’ 

corresponds  to  a  spurious  nonpropagating  mode  of  the  form  ■  Since  Ai  >  2 

for  all  2;  =  |A;|cAt/v^  >  0,  the  mode  rapidly  decays  and  poses  no  threat  to  stability. 

The  other  two  roots  are  a  pair  of  complex  conjugates: 


A2/3  =  -  ^(1  T  iV3)  (27 z^  +  3V3^/27z^  +  2^2  + 

1  ±*^3 


1/3 


6(27:22  +  3y3v/27FT^  +  ly/^ 


4/3 


l±tV2z-z^T  +  4:2"  ±  +  0(zy  as  z  ^  0. 


We  can  show  that 


2  16  ,  -  16W3  -  4W2  -  16W  +  4 


|A2/3|^  =  f  + 


36  W2 
T/3 


>  1  for  all  hh  >  1 


where  W  =  {27z‘^  +  3\p3\/27z'^  +  2^2  -|-  >1  for  all  z  >H. 

Since  A  =  we  have 

1 


UJ  = 


iAt 

1 


log(A) 


-  log(l  +  iV2z  -z^-  'TAz?  +  4^1  +  S^^z}  +  0(zZ)) 

« liic  (1  -  _  POAAA  +  0((|i|cAi)-^: 

So  the  phase  error  is  \-^  —  1|  =  0((\k\cAty). 


The  presence  of  the  imaginary  third  term  in  the  expansion  shows  that  ~  ^-in(-i{\k\cAtf)At) 
^-n{\k\c)^ At'^ /2 ^  causing  the  mode  to  decay.  This  is  why  we  term  this  scheme  diffusive  (or  dis¬ 
sipative). 
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7.2.2.  Dispersive  Scheme.  Substitute  ansatz  u'^  =  u)n/\t) 


+  2u’^  +  =  Au 


Obtain  a  polynomial 


A2  +  2 


\ \k\^  +  ) 


A  +  1 


0 


where  A  = 

We  can  solve  to  obtain  A1/2  =  ,  which  gives  IA1/2I  = 

non-dissipative. 

Noting  that  cos(AAt)  =  |(Ai  +  A2),  and  defining  2;  =  \k\/a  = 


1,  meaning  this  scheme  is 
|fc|cAt/2,  we  obtain 


bj  =  arccos 
At 


|fcp  + 


=  ^  arccos 


At 


1  +  ^2 


—  (z  -  zy3  +  0{z^)) 


cAt)^  +  0((|fc|cAt)^ 


So  the  phase  error  is  ||||^  —  1|  =  0{{\k\cAt)‘^).  Moreover,  A  is  real,  owing  to  the  non- 
dissipative  nature  of  the  scheme. 


7.3.  Fully  Discrete  Von  Neumann  Analysis.  In  this  section,  we  provide  fully  discrete 
Von  Neuman  analyses  for  the  two  fully  discrete  schemes  derived  in  3,  and  show  that  they  are 
unconditionally  stable,  in  the  sense  of  Von  Neumann  analysis.  Combining  the  quadrature 
rules  and  exponential  recursion,  and  ignoring  boundaries,  we  can  write 

~  ^h[fj]  =  Pfj  +  -Q  {fj+1  +  fj-i)  +  R  {fj+i  —  ‘2fj  +  fj-i)  + 

^  00 

+  2  ^  ifj+k  +  fj-k)  +  Q  U'j+k+l  +  fj-k-l)  + 

RR  ifj+k+l  —  ‘^fj+k  +  fj+k-1  +  fj-k+1  —  “2 fj-k  +  fj-k-l)] 

with  u  =  ah  =  aAx  and  P,  Q,  and  R  defined  as  in  Section  3. 

Using  the  discrete  convolution  operator  Ih,  and  ignoring  sources  and  boundaries,  we  can 
write  the  diffusive  version  of  the  fully  discrete  scheme  as 

u]^^  =  hk[5u]-4u]-^  +  u]-^] 

and  the  dispersive  scheme  as 

+  2u]  +  =  4h[u]] 
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Defining  Ih,x  and  Ih,y  as  the  discrete  convolution  operators  acting  in  the  x-  and  y-directions, 
and  again  ignoring  sources  and  boundaries,  we  can  write  the  diffusive  version  of  the  2D  fully 
discrete  scheme  as 

“"f  = 

and  the  2D  fully  dispersive  scheme  as 

=  44,,|4,,|«”J] 

We  will  refer  to  these  as  free-space  schemes.  We  now  state  a  stability  theorem  based  on 
the  Von  Neumann  analysis  of  the  schemes. 

Theorem.  The  fully  discrete  dispersive  and  diffusive  free-space  schemes  in  one  dimen¬ 
sion  are  unconditionally  stable.  In  higher  dimensions,  the  corresponding  dimensionally-split 
schemes  are  also  unconditionally  stable. 

To  prove  the  stability  theorem,  we  consider  some  properties  of  the  discrete  covolution 
operator. 

Lemma.  Define  the  amplification  factor  A(k,u)  =  e  Then  A{k,  u)  is  well- 

defined  (independent  off),  and  the  following  hold. 

•  If  u  >  0  and  0  <  |fcAa;|  <  tt,  then  0  <  A{k,  u)  <  1. 

•  If  n  >  0,  then  2l(0,  z/)  =  1. 

•  For  any  0  <  \kAx\  <  n,  limj,_,.o+  A{k,  u)  =  0. 

•  For  any  0  <  \kAx\  <  n,  lim,^_,.oo  A{k,  z/)  =  1. 

Proof  of  the  Lemma.  We  calculate: 


A{k,u) 


=p(l  +  ^e  cos{kkAx)  + 


k=l 

oo 


Q  e'^  ^  e  cos{kkAx)  + 


k=l 


2R{cos{kAx)  —  1)  j  1  +  cosfikkAx) 

k=l 


=1  +  T, 


rp  ^-^{cos{kAx)  —  1)^  —  2^{dcos(kAx)  —  l)(cos(fcAa;)  —  1) 

—  2  cos{kAx)d  +  1 

with  d  =  e~'^. 

If  kAx  =  0  and  z/  >  0,  then  T  =  0.  If  0  <  \kAx\  <  tt,  some  calculus  shows  — 1  <  T  <  0 
for  any  z/  >  0,  and  lim;^_,.o+  T  =  —  1,  and  hmi,_,.oo  T  =  0.  □ 
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Proof  of  the  Theorem.  We  first  prove  that  each  one-dimensional  free-space  fully  discrete 
scheme  is  unconditionally  stable,  then  describe  the  extension  of  the  result  to  multiple  di¬ 
mensions  for  the  dimensionally-split  schemes. 

Diffusive  Scheme 

Consider  applying  the  fully  discrete  diffusive  scheme  to  uf  =  obtain  a 

polynomial 

+  5A  -  =  0 

where  A  =  and  =  2/A{k,i').  Note  that  the  lemma  implies  z  >  2  for  all  k  and 

Ax,  At  >  0.  The  condition  for  stability  is  |A|  >  1  for  all  roots  of  the  polynomial,  which  we 
verify  below. 

The  hrst  root  corresponds  to  a  spurious  mode; 


Ai 


We  can  show  that 


4  -  104z  +  10b  +  27z  -  52 

3  ^  3^ 

3\/3V3V27z^  -  lOiz  +  l00  +  27 z  ^ 


Ai 


W^  +  4W  +  1 

3W 


>  2 


for  all  W  >  1 


where  W  =  >  1  for  e  >  2. 

V2 

The  other  roots  are  pair  of  complex  conjugates: 


,  4  ,  ^373^2^  -  104z  +  100  +  27z  -  52 

A./3=3-(1T.^) - - 


-  (1 


6^3v/3v^2  _  104;^  +  100  +  27;^  -  52 


We  can  show  that 

IA2/3P  = 


2  4W^  -  16W^  +  60W2  -  16W  +  4 


36W2 


>  1  for  all  W  >  1 


where  W  =  >  1  for  e  >  2. 

This  proves  the  theorem  in  the  case  of  the  one  dimensional  diffusive  scheme. 


Dispersive  Scheme 

In  this  section,  we  prove  that  the  one-dimensional  free-space  fully  discrete  dispersive  scheme 
is  unconditionally  stable  and  non-dissipative.  Consider  applying  the  fully  discrete  dispersive 
scheme  to  uf  =  ^  obtain  a  polynomial 

A^  +  (2  -  4:z)\  +  1  =  0 

where  A  =  and  =  A{k,u).  Note  that  the  lemma  implies  0  <  <  1  for  all  k  and 

Ax,  At  >  0.  We  can  solve  to  obtain  the  roots  A1/2  =  {2z  —  1)  ±  —  z)z.  We  can  show 

IA1/2I  =  1  for  all  0  <  <  1,  so  the  fully  discrete  dispersive  scheme  is  unconditionally  stable, 

and  non-dissipative.  This  proves  the  theorem  in  the  case  of  the  one  dimensional  dispersive 
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scheme. 


Extension  to  Higher  Dimensions 

When  applying  the  dimensionally-split  two-dimensional  schemes  to 
we  obtain  the  same  Von  Neumann  polynomials  as  in  the  ID  case,  and  basically  the  same  sta¬ 
bility  analysis  can  be  repeated.  This  can  be  easily  extended  to  dimensionally-split  schemes 
in  any  dimension. 

□ 

This  Von  Neumann  analysis  does  not  take  into  consideration  the  effects  of  boundary 
conditions,  and  in  principle,  certain  numerical  boundary  conditions  could  result  in  instability. 
In  the  test  problems  presented  in  [8,  7,  9]  as  well  as  the  present  work,  the  stability  of 
the  method  seems  robust  under  a  variety  of  numerical  boundary  conditions.  A  stability 
analysis  for  some  ID  fully  discrete  schemes  (slightly  different  from  those  presented  here) 
with  numerical  Dirichlet  and  Neumann  boundary  conditions  was  carried  out  in  [8],  showing 
unconditional  stability  in  those  schemes.  A  similar  stability  analysis  could  also  be  carried 
out  for  the  schemes  considered  in  this  work  to  study  the  stability  of  these  schemes  under  the 
inclusion  of  numerical  boundary  conditions. 

References 

1.  TD  Arber  and  RGL  Vann,  A  critical  comparison  of  eulerian-grid-hased  vlasov  solvers,  Journal  of  com¬ 
putational  physics  180  (2002),  no.  1,  339-357. 

2.  Kiyoshi  Asano  and  Seiji  Ukai,  On  the  vlasov-poisson  limit  of  the  vlasov-maxwell  equation,  Studies  in 
Mathematics  and  Its  Applications  18  (1986),  369-383. 

3.  Willard  H  Bennett,  Magnetieally  self-foeussing  streams.  Physical  Review  45  (1934),  no.  12,  890. 

4.  Charles  Kennedy  Birdsall  and  Allan  Bruce  Langdon,  Plasma  physies  via  eomputer  simulaition,  CRC 
Press,  2005. 

5.  Jose  A  Bittencourt,  Fundamentals  of  plasma  physics.  Springer,  2004. 

6.  Oscar  P  Bruno  and  Mark  Lyon,  High-order  unconditionally  stable  fc-ad  solvers  for  general  smooth  do¬ 
mains  i.  basic  elements.  Journal  of  Computational  Physics  229  (2010),  no.  6,  2009-2033. 

7.  M.  Causley,  Y.  Ciiglii,  E.  Wolf,  and  A.  Christlieb,  A  Fast,  Unconditionally  stable  solver  for  the  wave 
equation  based  on  the  Method  of  Lines  Transpose,  (in  preparation). 

8.  Matthew  Causley,  Andrew  Christlieb,  Benjamin  Ong,  and  Lee  Van  Croningen,  Method  of  lines  transpose: 
An  implieit  solution  to  the  wave  equation.  Mathematics  of  Computation  (2014). 

9.  Matthew  F  Causley  and  Andrew  J  Christlieb,  Higher  order  a-stable  sehemes  for  the  wave  equation  using 
a  suceessive  eonvolution  approaeh,  SIAM  Journal  on  Numerical  Analysis  52  (2014),  no.  1,  220-235. 

10.  Roman  Chapko,  Rainer  Kress,  et  ah,  Rothes  method  for  the  heat  equation  and  boundary  integral  equations. 
Journal  of  Integral  Equations  and  Applications  9  (1997),  no.  1,  47-69. 

11.  Cuangye  Chen,  Luis  Chacon,  and  Daniel  C  Barnes,  An  energy-  and  charge- conserving,  implicit,  electro¬ 
static  particle-in- cell  algorithm.  Journal  of  Computational  Physics  230  (2011),  no.  18,  7018-7036. 

12.  Chio-Zong  Cheng  and  Ceorg  Knorr,  The  integration  of  the  vlasov  equation  in  configuration  space.  Journal 
of  Computational  Physics  22  (1976),  no.  3,  330-351. 

13.  Andrew  J  Christlieb,  Robert  Krasny,  John  P  Verboncoeur,  Jerold  W  Emhofl,  and  Iain  D  Boyd,  Grid-free 
plasma  simulation  teehniques.  Plasma  Science,  IEEE  Transactions  on  34  (2006),  no.  2,  149-165. 

14.  Bruce  I  Cohen,  A  Bruce  Langdon,  and  Alex  Friedman,  Implicit  time  integration  for  plasma  simulation. 
Journal  of  Computational  Physics  46  (1982),  no.  1,  15-38. 

15.  Lawrence  C.  Evans,  Partial  differential  equations,  2nd  ed  ed.,  Graduate  studies  in  mathematics,  vol.  v. 
19,  American  Mathematical  Society,  Providence,  R.L,  2010. 

30 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


16.  Francis  Filbet  and  Eric  Sonnendriicker,  Comparison  of  eulerian  vlasov  solvers,  Computer  Physics  Com¬ 
munications  150  (2003),  no.  3,  247-266. 

17.  Matthew  R  Gibbons  and  Dennis  W  Hewett,  The  darwin  direct  implicit  particle-in- cell  (dadipic)  method 
for  simulation  of  low  frequency  plasma  phenomena,  Journal  of  Computational  Physics  120  (1995),  no.  2, 
231-247. 

18.  Roger  W  Hockney  and  James  W  Eastwood,  Computer  simulation  using  particles,  CRC  Press,  1988. 

19.  Mary  Catherine  A  Kropinski  and  Bryan  D  Quaife,  Fast  integral  equation  methods  for  rothe’s  method 
applied  to  the  isotropic  heat  equation.  Computers  &  Mathematics  with  Applications  61  (2011),  no.  9, 
2436-2446. 

20.  _ ,  Fast  integral  equation  methods  for  the  modified  helmholtz  equation.  Journal  of  Computational 

Physics  230  (2011),  no.  2,  425-434. 

21.  Lev  Davidovich  Landau  and  Evgenii  Mikhailovich  LifshitC  sCi,  The  classical  theory  of  fields,  vol.  2, 
Butterworth-Heinemann,  1975. 

22.  A  Bruce  Langdon,  Bruce  I  Cohen,  and  Alex  Friedman,  Direct  implicit  large  time- step  particle  simulation 
of  plasmas.  Journal  of  Computational  Physics  51  (1983),  no.  1,  107-138. 

23.  Jongwoo  Lee  and  Bengt  Fornberg,  Some  unconditionally  stable  time  stepping  methods  for  the  .3d 
maxwell’s  equations.  Journal  of  Computational  and  Applied  Mathematics  166  (2004),  no.  2,  497-523. 

24.  Peijun  Li,  Hans  Johnston,  and  Robert  Krasny,  A  cartesian  treecode  for  screened  coulomb  interactions. 
Journal  of  Computational  Physics  228  (2009),  no.  10,  3858-3868. 

25.  Mark  Lyon  and  Oscar  P  Bruno,  Fligh-order  unconditionally  stable  fc-ad  solvers  for  general  smooth  do¬ 
mains  ii.  elliptic,  parabolic  and  hyperbolic  pdes;  theoretical  considerations.  Journal  of  Computational 
Physics  229  (2010),  no.  9,  3358-3381. 

26.  PJ  Mardahl  and  JP  Verboncoeur,  Charge  conservation  in  electromagnetic  pic  codes;  spectral  comparison 
of  boris/dadi  and  langdon-marder  methods.  Computer  physics  communications  106  (1997),  no.  3,  219- 
229. 

27.  Martin  Masek  and  Paul  Gibbon,  Mesh-free  magnetoinductive  plasma  model.  Plasma  Science,  IEEE 
Transactions  on  38  (2010),  no.  9,  2377-2382. 

28.  Takefumi  Namiki,  3-d  adi-fdtd  method-unconditionally  stable  time-domain  algorithm  for  solving  full  vec¬ 
tor  maxwell’s  equations.  Microwave  Theory  and  Techniques,  IEEE  Transactions  on  48  (2000),  no.  10, 
1743-1748. 

29.  Donald  W  Peaceman  and  Henry  H  Rachford,  Jr,  The  numerical  solution  of  parabolic  and  elliptic  differ¬ 
ential  equations.  Journal  of  the  Society  for  Industrial  &  Applied  Mathematics  3  (1955),  no.  1,  28-41. 

30.  SV  Petropavlovsky  and  SV  Tsynkov,  Quasi-lacunae  of  maxwell’s  equations,  SIAM  Journal  on  Applied 
Mathematics  71  (2011),  no.  4,  1109-1122. 

31.  E  Pohn,  Magdi  Shoucri,  and  G  Kamelander,  Eulerian  vlasov  codes.  Computer  physics  communications 
166  (2005),  no.  2,  81-93. 

32.  RJ  Procassini,  CK  Birdsall,  and  EC  Morse,  A  fully  kinetic,  self-consistent  particle  simulation  model  of 
the  collisionless  plasma-sheath  region.  Physics  of  Fluids  B:  Plasma  Physics  (1989-1993)  2  (1990),  no.  12, 
3191-3205. 

33.  Erich  Rothe,  Zweidimensionale  parabolische  randwertaufgaben  als  grenzfall  eindimensionaler  randwer- 
taufgaben,  Mathematische  Annalen  102  (1930),  no.  1,  650-670. 

34.  Jack  Schaeffer,  The  classical  limit  of  the  relativistic  vlasov-maxwell  system.  Communications  in  mathe¬ 
matical  physics  104  (1986),  no.  3,  403-421. 

35.  David  N  Smithe,  John  R  Cary,  and  Johan  A  Carlsson,  Divergence  preservation  in  the  adi  algorithms  for 
electromagnetics.  Journal  of  Computational  Physics  228  (2009),  no.  19,  7289-7299. 

36.  Alen  Taflove  and  Susan  C  Hagness,  Computational  electrodynamics:  the  fdtd  method,  Artech  House 
Boston,  London  (2000). 

37.  John  P  Verboncoeur,  Aliasing  of  electromagnetic  fields  in  stair  step  boundaries.  Computer  physics  com¬ 
munications  164  (2004),  no.  1,  344-352. 

38.  _ ,  Particle  simulation  of  plasmas:  review  and  advances.  Plasma  Physics  and  Controlled  Fusion  47 

(2005),  no.  5A,  A231. 

39.  Kane  S  Yee,  Numerical  solution  of  initial  boundary  value  problems  involving  maxwell’s  equations,  IEEE 
Trans.  Antennas  Propag  14  (1966),  no.  3,  302-307. 

31 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


40.  Fenghua  Zhen,  Zhizhang  Chen,  and  Jiazong  Zhang,  Toward  the  development  of  a  three-dimensional 
unconditionally  stable  finite- difference  time-domain  method,  Microwave  Theory  and  Techniques,  IEEE 
Transactions  on  48  (2000),  no.  9,  1550-1558. 

41.  Fenghua  Zheng,  Zhizhang  Chen,  and  Jiazong  Zhang,  A  finite- difference  time-domain  method  without  the 
courant  stability  conditions,  Microwave  and  Guided  Wave  Letters,  IEEE  9  (1999),  no.  11,  441-443. 

Department  of  Mathematics,  Michigan  State  University,  East  Lansing,  MI  48824 
E-mail  address:  christli@msu.edu 


32 


DISTRIBUTION  A:  Distribution  approved  for  public  release. 


Time  Step  n  =  1 


Time  Step  n  =  1 


Particle  Position  x 


'o 

o 

Q> 

> 

a) 

o 

tr 

C3 

CL 


Particle  Position  x 


(A) 


(b) 


Time  Step  n  =  500 


>s 

2 

'o 

o 

0 

> 

0 

0 

o 

r 

-2 

0 

Q_ 

5 

Particle  Position  x 

(e) 


Time  Step  n  =  1000 


Particle  Position  x 

(G) 


Time  Step  n  =  500 


'o 

_o 

0 

> 

o 

tr 

C3 

CL 


-5 


Particle  Position  x 

(f) 


Time  Step  n  =  1000 


'o 

_o 

0 

> 

_0 

o 

c5 

CL 


Particle  Position  x 

(H) 


Figure  5.  We  see  selected  particle  phase  space  plots  for  the  two  stream  instability 
problem.  The  left  column  is  from  the  ID  simulation,  while  the  right  column  is 
from  the  2D  simulation,  following  a  fixed  slice  of  particles  initialized  along  the  line 

y  — 
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Figure  6.  Landau  damping  of  the  lowest  mode,  corresponding  to  k  =  0.5. 
Green  is  the  ID  numerical  result,  blue  is  the  2D  numerical  result  (measured 
along  the  central  y  =  0  slice),  and  red  is  the  prediction  of  linear  theory.  We 
see  that  the  correct  decay  rate  is  reproduced  in  our  simulations 
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Figure  7.  In  7a,  we  see  the  scalar  potential  profile  at  t  =  3.6  thermal-ion  transit 
times.  In  7b,  we  see  the  simulation  electron  and  ion  count,  the  red  and  blue  curves 
respectively,  along  with  the  injection  rate,  the  black  dashed  line. 
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.  The  staggered  grid  used  for  the  Bennett  pinch  problem. 
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Figure  9.  The  figure  shows  9a  electron  density,  9b  and  potential  energies,  9c 
magnetic  field  at  various  CFL  numbers,  and  9d  the  relative  error  in  the  azimuthal 
magnetic  field  Bq  (normalized  to  the  maximum  value  of  the  magnetic  field) .  Results 
are  with  a  CFL  number  of  3,  except  as  noted  in  9c.  Position  units  are  in  terms  of 
effective  beam  radius  Rh- 
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Figure  10.  The  staggered  grid  used  for  the  Mardahl  beam  problem. 
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Figure  11.  The  figure  shows  the  divergence  error  in  the  electric  fields  and  the 
final  beam  distribution  calculated  from  a  wave  equation  potential  11a,  11c  and  in 
the  Poisson  equation  potential  lib,  lid. 
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